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ABSTRACT 

An  efficient  longwave  radiation  transfer  band  model  was 
developed  and  applied  to  the  MIT  Stratospheric  General  Cir¬ 
culation  Model  (GCM) .  The  band  model  is  a  revised  version 
of  models  by  Ramanathan  (1976)  and  by  Alimandi  and  Visconti 
(1979) .  It  is  validated  by  comparison  with  work  by  Dopplick 
(1971) ,  Slade  (1969)  and  by  Park  and  London  (1974)  . 

The  band  model  produced  pole-to-pole  heating  gradients 
up  to  9  degK/day  compared  with  Newtonian  heating  gradients 
of  1-2  degK/day  used  previously  in  the  GCM.  In  5  day  model 
integrations,  the  band  model  sharpened  erroneous  tempera¬ 
ture  gradients  in  the  GCM  when  compared  with  the  Newtonian 
heating  scheme.  Poleward  eddy  heat  fluxes  appear  to  be  too 
weak  in  the  low-resolution,  6  wave  spectral  GCM  to  balance 
the  heating  gradients  of  the  band  model.  An  18-wave  spectral 
GCM  is  planned  by  the  MIT  stratospheric  working  group. 

With  the  band  model,  global  heating  rates  up  to  3  degK/ 
day  are  noted.  Such  global  heating  rates  cannot  be  sustained 
in  the  atmosphere.  Methods  to  achieve  global  heating  balance 
are  explored.  The  most  straight  forward  method  is  to  set  the 
horizontal-averaged  temperatures  of  the  GCM  equal  to  the 
radiative  equilibrium  temperatures  of  the  band  model. 


Thesis  Advisor:  Dr.  Ronald  G.  Prinn,  Associate  Professor  of 

Meteorology  and  Physical  Oceanography 
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1.  INTRODUCTION 

A  numerical  simulation  of  the  general  circulation  must 
assume  some  knowledge  of  the  diabatic  sources  and  sinks  of 
heat  energy.  How  much  detail  of  the  radiative  energy  dis¬ 
tribution  can  be  included  in  a  particular  general  circu¬ 
lation  model  (GCM)  depends  on  the  problem  considered,  the 
accuracy  desired  and  the  runtime  constraints  required. 

This  paper  addressed  the  need  for  efficient  radiative 
transfer  computations  in  the  stratosphere  and  mesosphere 
when  accurate  simulation  is  sought  under  runtime  constraints. 

Murgatroyd  and  Goody  (1958)  presented  one  of  the  first 
summaries  of  net  diabatic  heating  rates  due  to  radiation 
in  the  stratosphere  and  mesosphere.  Dopplick  (1971) ,  Slade 
(1969) ,  and  Tahnk  and  Newell  (1975)  have  reported  heating 
rates  due  to  the  major  stratospheric  absorbing  gases, 

CC>2  and  ,  using  automated  radiation  band  models.  Park 
and  London  (1974)  include  the  minor  stratospheric  gas,  H20. 
Dickinson  (1972)  made  revisions  to  the  CC>2  band  model  by 
including  the  effects  of  minor  isotopic  and  "hot"  bands. 

Leovy  (1964)  developed  theoretical  dynamical  models 
for  the  upper  atmosphere  based  on  an  idealized,  symmetrical 
version  of  the  radiative  heating  distribution.  Manabe  and 
Hunt  (1968)  incorporated  radiative  band  models  into  a  stra¬ 
tospheric  general  circulation  model  up  to  40  km.  Dickin¬ 
son  (1973)  addressed  the  need  for  rapid  radiative  computa¬ 
tions  for  use  in  GCMs,  suggesting  a  Newtonian  cool-to-space 
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approximation  for  infrared  cooling.  This  method  takes 
negligible  computation  time  but  neglects  layer-to- layer 
exchanges  of  heat. 

Ramanathan  (1976)  reported  a  radiative  transfer 
band  model  which  he  claimed  to  be  an  order  of  magnitude 
faster  than  the  detailed  radiative  models  by  Dickinson 
(1972)  and  Ellingson  (1972)  and  yet  accurate  to  within  1% 
of  their  results.  In  1978-79  Visconti  (1979)  modified  the 
long  wave  radiation  portions  of  the  Ramanathan  model  for 
the  gases  CC>2  and  C>3 .  The  purpose  of  the  work  reported 
here  is  to  further  streamline  this  radiative  transfer  code 
and  incorporate  it  into  the  MIT  stratospheric  general  cir¬ 
culation  model  where  it  replaces  a  Newtonian  cooling-to- 
space  scheme. 

To  date,  the  computation  time  for  CC>2  has  been  re¬ 
duced  by  a  factor  of  20 .  The  ozone  computation  time  has 
been  reduced  by  a  factor  of  two.  These  two  components  have 
been  applied  to  the  GCM  above  30  km.  Net  radiative  heating 
rates  compare  well  with  Park  and  London  (1974) .  Increases 
in  magnitude  up  to  a  factor  of  6  in  the  pole-to-pole 
heating  gradient  are  found  in  comparison  with  the  Newtonian 
cooling  approximation. 

The  MIT  GCM  predicts  a  temperature  gradient  directed 
from  the  winter  pole  to  the  summer  pole  at  70  km  in  contrast 
with  observations.  The  initial  effect  of  the  more  realistic 
radiative  heating  rates  is  to  sharpen  the  winter-summer 
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gradient  in  the  model. 

Horizontally-averaged  heating  rates  up  to  +3  deg/ 
day  (30%  of  the  maximum  pole  to  pole  differential  heating) 
are  found  at  upper  atmospheric  levels.  This  net  heat 
imbalance  appears  because  the  band  model  radiative  equili¬ 
brium  temperature  profile  is  warmer  than  the  reference 
temperature  profile  of  the  GCM.  An  empirical  parametriza- 
tion  of  sub-grid  scale  vertical  eddy  heat  diffusion  is 
presented  as  one  method  to  accommodate  the  heat  imbalance. 


2.  RADIATION  THEORY  AND  THE  RAMANATHAN/VISCONTI  BAND 

MODEL 

In  the  long  wave  range  of  the  electromagnetic  spectrum, 
radiant  energy  is  produced  by  emission  everywhere  in  the 
atmosphere.  For  a  given  frequency,  v,  the  flux  of  radiant 
energy  at  a  point  is  the  integral  of  the  emitted  energy 
from  every  point  in  the  atmosphere  attenuated  according  to 
Beer -Lambert 's  Law  of  exponential  decay  through  absorption 
by  the  intervening  absorbing  gases.  Long  wave  cooling,  Q, 
is  the  result  of  a  divergence  of  this  three  dimensional 
flux.  (See,  for  more  detail.  Goody,  1964) 


Qv 


(i) 


p  =  density 

0^  -  constant  pressure  heat  capacity 
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Atmospheric  band  models  generally  simplify  the 
geometry  of  radiative  transfer  by  assuming  a  vertically 
stratified  plane  parallel  atmosphere  with  horizontally 
homogeneous  temperature  and  absorber  gas  compositions 
at  each  grid  point  considered.  Integration  over  zenith 
angle  can  be  parameterized  in  the  plane  parallel  atmosphere 
with  a  constant  "dif fusivity "  factor  so  that  the  problem 
is  reduced  to  one  dimension: 


Q 


v 


(2) 


Rodgers  and  Walshaw  (1966)  compiled  evidence  that  the 
cooling  rate  error  due  to  using  a  constant  dif fusivity  is 
1.5  per  cent  or  less  when  compared  with  more  complex 
transmission  dependent  zenith-integrated  dif fusivities . 

In  the  Ramanathan  model,  the  carbon  dioxide  15y  band 
and  the  ozone  9.6y  band  are  each  handled  by  empirically 
derived  band  absorptance  formulations ,  in  which  a  mean 
absorptance  for  an  entire  band  is  defined  for  a  given 
atmospheric  path  based  on  temperatures,  pressures,  and 
absorbing  gas  amounts  along  the  path  together  with  other 
predetermined  parameters. 

Following  Manabe  and  Moller  (1961)  ,  the  appropriate 
version  of  the  radiative  transfer  equation  for  F(z) ,  the 
net  flux  at  a  given  level,  is  obtained  by  integrating  the 
basic  equation  (3)  vertically  from  the  ground  to  level  Z 
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and  then  from  the  top  of  the  atmosphere  to  level  Z  with 
black  body  radiation  at  the  ground  and  no  downward 
radiation  at  the  top  of  the  atmosphere . 


dl (Z') 

dx  =  I(Z')  -  J(Z')  (3) 


where  I  is  the  radiation  intensity  along  a  path 
J  is  the  source  function  for  emission  and 
x  is  the  optical  depth 

Using  the  Planck  function,  B,  as  the  source  function,  multi- 

—  T 

plying  by  the  integrating  factor, T  =  e  ,  and  integrating 
by  parts  gives  the  upward  flux,  I  +  (Z)  and  downward  flux, 

I+(Z)  : 

Zsf  c 

I .  ( Z)  =  B  { Z)  +  /  r  ( Z  ' - S )  dB 

f  z 

I+(Z)  =  B(Z>Z  -  B(ztop)r(ztop-z') 

+  /  t0pr(Z'-Z)dB  (4) 

z 

Introducing  the  absorptance,  A  =  1  -  r,  the  net  flux  at 
level  Z  is  given  by  equation  (5) 


Net  Flux,  F(Z)  =  I+  -  I+  =  B(Zsfc)  -  B ( ZtQp) A(Ztop, Z) 


^sf  c 

-  /  A ( Z ' , Z) dB ( Z  ' ) 
Ztop 


(5) 


where  A(Z',Z)  represents  the  absorptance  between  level 

Z'  and  level  Z. 

The  above  expression  is  exact  when  applied  to  a  single 
wave  frequency.  V7hen  equation  (5)  is  applied  to  a  band, 
the  absorptance  must  be  parameterized  empirically.  Equations 
(2)  and  (5)  are  the  basic  cooling  rate  and  flux  computa¬ 
tion  equations  used  in  the  Ramanathan  model. 

The  15yCC>2  and  9.6yC>3  bands  are  sufficiently  narrow 
that  a  single  wave  frequency  can  be  used  to  evaluate  the 
Planck  function:  667  cm  ^  for  CC>2  and  1041  cm"^  for  . 

The  empirical  absorptance  formulations  developed  by 
Cess  and  Tiwari  (1972)  are  used  by  Ramanathan.  For  carbon 
dioxide,  the  pressure-broadened  absorptance  is  given  by 

io  li 

A(£)  =  2A  Jtn{  1  +  Z  c.  }  (6) 

°  i=l  1 

4v 

where  =  1.66  (v— /  S.q.P  PdZ 
l  AD.  lie 

o  l 

VQ  =  mean  line  half  width  spacing 
Aq  =  band  width  parameters 
q^  =  relative  abundance  of  isotopes  to  the 
primary  isotope 
=  mean  line  spacing 
Si  =  band  intensity 

Pc  =  partial  pressure  of  carbon  dioxide 
P  =  atmospheric  pressure 
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where  four  isotopic  and  six  "hot"  bands  are  considered. 

A  strong  line  approximation  has  been  used. 

Alimandi  and  Visconti  (1979)  achieved  a  reduction 
in  computation  time  by  noting  that  the  band  width  param- 

-1  h 

eter,  Aq,  may  be  scaled  as  Aq  =  22.3  cm  (T/300)  so  that 
A^ j ,  the  absorptance  between  two  pressure  levels  i  and  j, 
can  also  be  scaled  like  the  square  root  of  temperature  at 
the  lower  level,  i. 


A. .  =  A. .ref (T./T.ref) * 
l  l 


(7) 


The  band  absorptance  formula  for  Doppler  broadening 
used  by  Ramanathan  is  shown  in  equation  (8) . 


A(u)  =  Aqu(1  -  0.18U/S)  u/6 <1.5  (8) 

=  0.753  Aq6{  [ln(u/6)  )3/2+1.2]}  u/6  >  1.5 


where  <5  =  Kvd 

=  Doppler  line  half  width 
u  =  optical  path  length 


=  —  P 
Ao  c 

Both  the  doppler  half  width,  v^,  and  the  band  width 
parameter,  Ao»  are  proportional  to  the  square  root  of 
temperature  so  for  u/6  >  1.5,  doppler  absorptance  scales 
like  temperature: 
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A.  .ref (T . /T . ref v 
13  a/  l  ) 


(9) 


Alimandi  and  Visconti's  version  of  the  CC>2  band  model 
computes  Doppler  and  Lorentz  absorptances  for  a  reference 
temperature  profile  and  uses  the  scaling  described  to 
compute  absorptances  for  each  new  temperature  profile. 

Using  a  seasonal  and  latitudinal  range  of  temperature  pro¬ 
files,  Alimandi  and  Visconti  (1979)  found  accuracy  remained 
within  10%  of  the  exact  computation  when  the  scaling  was 
used . 

For  ozone,  the  strong  line  approximation  is  not  valid 
in  the  atmosphere.  Therefore  a  more  complex  band  absorp- 
tance  formulation  must  be  used  for  pressure  broadening. 

Also,  ozone  mixing  ratios  vary  with  season  and,  above  60  km, 
even  diumally.  For  these  reasons,  scaling  was  not  feasible 
in  the  03  band  model  and  the  absorptances  were  calculated 
for  each  profile  using  equaton  (10) . 

A  (u,  B)  =  2A  Jin  {1  +  - - - 1~}  (10) 

°  [4+u(i+i/e)r 

4v 

where  the  line  width  parameter,  (5  -  /  Pdu 

and  all  other  symbols  have  been  defined  in 
equations  (6)  and  (8). 
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3.  CODE  EFFICIENCY  IMPROVEMENTS 

Radiative  transfer  is  inherently  a  time  consuming 
computation.  To  take  account  of  exchange  of  radiative 
energy  between  the  layers  of  the  atmosphere  in  an  n  -  layer 
model  requires  that  code  be  repeated  at  least  n  x  n  times. 

In  fact  greater  vertical  resolution  may  be  required.  To 
use  the  Ramanathan  band  absorptance  formulation  required 
a  vertical  resolution  of  better  than  3  km.  For  the  MIT 
26  level,  0-7 2km  stratospheric  GCM  that  meant  the  flux  compu¬ 
tation  code  was  repeated  2  x  n  x  n  (1352)  times  for  a 
single  profile. 

The  second  feature  of  radiative  transfer  code  which 
makes  it  time  consuming  is  the  necessity  of  using  mathemati¬ 
cal  functions  such  as  the  exponential,  logarithm  and  square 
root  functions.  The  exponential  is  required  for  the  Planck 
function  of  a  single  frequency.  In  Ramanathan' s  band  ab¬ 
sorptance  formulation  both  natural  logarithm  and  square 
root  are  required. 

Efficiency  improvements  are  made  by  concentrating 
on  the  features  of  repeated  code  and  mathematical  functions. 
Table  1  lists  the  results  of  improvements  made  in  average 
run  time  for  one  time  step  on  December  1  in  the  MIT 
Stratospheric  GCM.  The  band  model  heating  code  is  used  at 
240  grid  points  for  incorporation  into  the  six-wave, 
spectral  GCM. 

For  the  CO 2  code  the  first  improvement  was  to  compute 
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and  table  the  Planck  function  B  =  and  scaling 

j,  e 

factors  T^/T^ref,[T^/T^ref  ]  ^  for  each  layer  of  the  band 
model.  Computation  time  dropped  by  a  factor  of  six  with 
no  loss  of  accuracy.  (See  Table  1,  Section  1.) 

The  second  improvement  relied  on  a  feature  of  the 
CDC  computer  hardware.  In  the  CDC  central  processor,  the 
instruction  word  stack  (IWS)  holds  up  to  48  instructions. 

If  a  do  loop  will  fit  completely  in  the  IWS,  computation 
is  much  more  efficient  since  instructions  need  not  be 
fetched  from  memory  in  each  repetition  of  the  loop.  For 
the  CO ^  flux  computation,  separating  the  upward  and  down¬ 
ward  flux  components  resulted  in  two  do  loops  each  of  which 
could  fit  in  the  IWS.  Because  the  flux  computation  domi¬ 
nates  the  code  time,  an  overall  speed-up  factor  of  two  was 
achieved,  still  with  no  loss  of  accuracy. 

Finally,  tables  of  the  Planck  function  and  square 
root  were  prepared  as  functions  of  temperature  in  the 
range  of  150-350°  K.  The  tables  were  computed  for  every 
degree  Kelvin.  The  square  root  table  gave  values  with  .1% 
accuracy  without  interpolation.  The  Planck  function  was 
accurate  within  1%  using  linear  interpolation.  A  further 
reduction  in  speed  of  30%  was  achieved.  Heating  rates 
remained  accurate  within  .1  deg/day. 

For  the  03  code,  precomputing  the  Planck  and  square 
root  functions  of  temperature  resulted  in  speeding  the  code 
by  a  factor  of  2.  (See  Table  1,  Section  2.)  This  was  not 


as  large  as  the  effect  in  the  CO 2  code.  The  primary 
reason  it  was  not  as  effective  here  was  because  the  band 
absorptance  calculation  still  required  a  logarithm  and 
square  root  to  be  performed  at  each  level  of  the  flux 
computation.  If  tables  of  these  functions  could  give 
sufficient  accuracy  another  factor  of  two  could  be  realized. 

Even  if  tables  were  successful  the  total  GCM  runtime 
would  be  doubled  by  the  code  alone.  Since  ozone  9.6y 
cooling  is  never  greater  than  20%  of  the  C02  cooling  and 
is  not  subject  to  rapid  variation,  less  frequent  calcula¬ 
tion  could  be  justified.  In  the  run  described  in  later 
sections  ozone  cooling  was  calculated  once  a  day  and  its 
effect  was  smoothed  over  24  hours  to  prevent  discontinui¬ 
ties.  This  was  effectively  a  drop  in  runtime  by  a  factor 
of  24. 

The  unaltered  radiation  code  overwhelmed  the  GCM  run¬ 
time  by  a  factor  of  10 .  The  described  revisions  fit  the 
infrared  radiation  code  into  30%  of  the  total  runtime. 

The  dominance  of  the  flux  computation  in  the  band 
model  was  clear.  In  the  fastest  C02  code,  the  flux  compu¬ 
tation  involved  only  five  multiplications,  a  conditional 
branch,  and  three  additions  for  each  layer-exchange. 

Still,  the  flux  computation  took  92%  of  the  C02  band  model 
time  (.6  sec/. 65  sec  per  time  step).  Thus  speeding  other 
portions  of  the  code  by  more  efficient  indexing,  identify¬ 
ing  invariant  code,  and  other  optimizations  did  not  lead  to 


significant  improvements  overall.  Further  major  improve¬ 
ments  will  require  a  reduction  in  the  number  of  layer  ex¬ 
changes  considered  at  each  time  step  in  the  flux  computa¬ 
tions  . 

A  sensitivity  analysis  was  performed  for  an  isolated 
5°  temperature  rise  at  each  level  of  the  model.  Results 
are  summarized  in  Figure  5  for  typical  levels.  Cooling 
rate  changes  of  more  than  .1  deg/day  occurred  at  the  level 
perturbed  and  at  the  next  higher  level  from  the  ground 
to  27  km.  Above  27  km,  .1  deg/day  changes  or  greater  were 
found  at  not  more  than  six  levels:  one  level  below  to 
four  levels  above  the  perturbed  level.  Thus  layer  to  layer 
exchange  need  only  be  computed  at  six  layers  on  an  hour  to 
hour  basis .  Flux  contributions  from  other  levels  can  be 
calculated  less  frequently  and  stored .  Another  factor  of 
5  in  computation  speed  would  be  saved  by  dropping  from  52 
half  levels  to  12  half-levels  contributing  to  calculations, 
but  a  baseline  of  the  time-integrated  affects  of  the  cur¬ 
rent  revision  of  the  model  should  be  set  first. 

4.  VALIDATION  OF  THE  MODEL 

For  validation  of  the  revised  radiative  band  model, 

I  used  ozone  and  temperature  analyses  given  by  Park  and 
London  (1974) .  Since  their  temperature  profile  does  not 
extend  to  the  surface ,  I  also  used  temperatures  from  Kuhn 
and  London  (Fig.  3,  1969).  For  near  surface  ozone  amounts. 
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I  used  values  from  an  MIT  GCM  run.  Park  and  London  do  not 
report  C02  and  cooling  separately  so  to  examine  these 
components  I  have  compared  them  with  results  from  Dopplick 
(1971)  and  Slade  (1969)  .  In  Figure  1/  15yC02  cooling  rates 
are  given.  In  agreement  with  Dopplick  and  Slade  the 
analysis  shows  maximum  cooling  at  the  stratopause  in  the 
summer  data.  Warming  occurs  near  the  equatorial  tropopause 
and  in  the  summer  data  at  70  km.  One  difference  is  that 
the  revised  Ramanathan  model  also  finds  warming  throughout 
the  data  below  3  to  8  km.  When  applied  in  the  MIT  GCM 
the  Ramanathan  model  will  not  be  used  in  the  troposphere. 

In  Figure  2,  9.6yC>3  cooling  rates  are  displayed. 
Maximum  cooling  is  centered  on  the  stratopause  and  extends 
to  all  but  the  winter  polar  night  latitudes.  Warming  is 
found  throughout  the  lower  stratosphere.  In  both  regions 
the  Ramanathan  model  gives  weaker  effects  than  the  Dopplick 
and  Slade  models  up  to  1  degK/day.  Considering  differences 
in  the  temperature  profiles  used,  agreement  is  fairly  good. 

Combining  the  C02  and  0^  band  model  with  the  MIT  GCM 
UV  heating  code  resulted  in  the  net  heating  rates  shown 
in  Figure  3 .  July  data  is  presented  on  the  left  of  the 
figure  to  match  the  published  data  by  Park  and  London  (1974). 
IN  the  summer  hemisphere  on  both  analyses,  heating  peaks 
at  10  deg/day  at  50  km  in  the  polar  region.  A  ridge  of 
strong  heating  extends  towards  the  equator  sloping  downwards 
to  40  km.  This  secondary  maximum  is  stronger  in  results 


based  on  Ramanathan 1 s  code. 

In  the  winter  hemisphere,  maximum  cooling  is  at 
60  km.  Park  and  London  have  a  higher  maximum  cooling  (-10 
deg/day  vs.  -7  deg/day) .  The  major  difference  in  the  anal¬ 
ysis  structure  appears  to  be  the  extension  of  the  cooling 
rate  "trough"  at  60  km  southward  across  the  analysis  in 
Figure  3.  Park  and  London's  data  do  not  show  a  feature  like 
this . 

Figure  4  is  an  analysis  of  the  GCM  UV  heating  rates. 
Comparison  wit*  Park  and  London's  UV  heating  rates  shows 
that  the  equatorial  cold  trough  and  enhanced  65-70  km 
heating  xn  the  summer  hemisphere  are  due  to  differences 
in  the  tv  heating  rates.  These  in  turn  are  probably 
because  the  published  ozone  analysis  is  for  daytime  amounts 
whi^e  the  net  UV  rates  published  were  integrated  around 
the  24  hour  cycle  during  which  ozone  amount  can  vary  sub¬ 
stantially  especially  at  high  altitudes. 

The  greater  cooling  in  the  winter  hemisphere  seen 
by  Park  and  London  is  due  at  least  in  part  to  the  inclu¬ 
sion  of  cooling  by  the  80  y  1^0  band.  Kuhn  and  London 
(1969,  Fig.  8)  show  radiative  cooling  rates  up  to  -1  deg/ 
day  for  a  water  vapor  mixing  ratio  of  10  6  in  the  upper 
atmosphere.  They  consider  10~6  the  most  probable  value. 
Estimating  numerical  values  from  the  published  temperature 
and  ozone  analyses  is  another  source  for  differences. 


5. 


DIABATIC  HEATING  IN  THE  MIT  STRATOSPHERIC  GCM 


The  MIT  Stratospheric  General  Circulational  Model 
was  designed  by  Cunnold  et  al.  (1975)  to  simulate  ozone 
production  and  transport.  Since  seasonal  ozone  amounts 
are  irregular,  several  years  of . integration  are  required 
to  reach  a  statistically  steady  state.  This  requirement 
led  to  the  choice  of  a  dynamic  "balance"  equation  system 
(Lorenz,  1960)  which  can  accept  a  time  step  up  to  an  hour. 

The  thermodynamic  prediction  equation  for  this 
system  is  given  in  equation  (11)  in  the  form  used  in  the 
GCM: 
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q'  =  diabatic  heating  per  unit  mass  minus  its 
horizontal  global  average 


In  this  equation,  only  deviations  from  the  reference 
temperature  profile  are  predicted.  To  assure  conservation 
of  energy,  average  temperature  on  a  model  level  is  not 
allowed  to  vary.  On  each  level,  the  static  stability 


parameter  is  held  constant.  Any  net  imbalance  in  the 
global  diabatic  heating  rate  for  a  level  is  removed 
equally  at  all  gridpoints .  These  features  can  potentially 
distort  the  radiative-dynamical  coupling.  Cunnold  et  al. 
(1975)  justify  these  approximations  as  not  affecting 
significant  conclusions  about  ozone  while  allowing  great 
improvements  in  computational  efficiency. 

Diabatic  heating,  q,  is  determined  empirically  in 
the  troposphere  so  as  to  reproduce  known  seasonal  estimates 
of  q  using  observed  seasonal  distributions  of  temperature. 
No  attempt  was  made  in  this  work  to  alter  this  portion 
of  the  heating  code. 

In  the  upper  atmosphere,  absorption  of  solar  radia¬ 
tion  by  ozone  is  computed  explicitly  as  the  product  of  the 
solar  constant,  1^ ,  times  the  absorption  coefficient  of 
ozone,  a  ,  ozone  mixing  ratio  and  the  transmission  func- 
tion,  e  T^,  all  integrated  over  wavelength  in  the  solar 
spectrum.  Infrared  cooling  is  parameterized  with  a  first 
term  Taylor  approximation  about  the  15pCC>2  cooling  rates 
for  the  1962  temperature  reference  profile.  The  cooling 
coefficients  used  are  due  to  Dickinson  (1973) .  This 
approximation  is  represented  as  Newtonian  cooling-to-space 
because  it  depends  only  on  the  local  temperature.  Radia- 
tional  exchange  with  other  levels  only  appears  in  the 
exact  cooling  rates  of  the  reference  profile.  Because 
this  scheme  is  linear  in  temperature,  the  calculation  could 
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be  done  in  spectral  coordinates  very  rapidly.  However, 
it  is  known  to  be  inexact  up  to  a  factor  of  2  at  70  km 
even  when  non-linear  pressure  and  temperature  corrections 
are  included.  (See  Dickinson,  1973,  Figure  4.) 

The  streamlined  Ramanathan  band  model  was  introduced 
into  the  model  in  grid  format  immediately  following  the 
UV  heating  calculation  which  is  also  done  in  grid  format. 
Conversion  to  spectral  form  for  use  in  the  prediction 
equations  is  done  using  a  fast  fourier  transform.  Spectral 
representation  leads  to  errors  up  to  1  degK/day  especially 
for  the  UV  heating.  UV  heating  rates  resemble  a  step 
function  at  the  edge  of  the  polar  night  before  spectral 
smoothing  and  are  therefore  hard  to  represent. 

Figure  6  displays  an  analysis  of  model  temperatures 
for  December  1  at  the  beginning  of  the  second  year  of  RUN 
34  in  the  MIT  stratospheric  GCM.  In  comparison  with  the 
analysis  of  observed  temperature  of  Park  and  London  (1974)  , 
the  model  analysis  shows  temperatures  which  are  too  cold 
throughout  the  polar  region  above  20  km.  This  feature 
has  been  typical  of  most  GCMs .  Figures  7  and  8  show  the 
net  heating  rates  on  December  1  using  respectively,  Newtonian 
cooling  and  the  revised  Ramanathan  band  model.  The  global 
mean  heating  has  been  removed  at  each  level  so  these 
figures  show  the  values,  q'/c  ,  used  in  the  prediction 

Sr 

equation  (11).  In  Figure  7  cooling  rates  are  generally 
small.  Newtonian  cooling  appears  to  match  closely  the  UV 


heating  everywhere.  Heating  is  strongest  in  the  equatorial 
stratopause  region.  Cooling  is  strongest  in  the  polar 
night  except  at  40  km  where  warming  occurs  at  the  temper¬ 
ature  minimum.  By  contrast.  Figure  8  shows  an  eight 
degree  pole-to-pole  gradient  at  55  km  primarily  due  to 
stronger  radiative  cooling  in  the  polar  night. 

Figures  10  and  11  show  the  effects  of  5  days  model 
integration  on  the  model  temperature  field  by  the  two 
cooling  schemes.  Comparing  Figures  6  and  10  the  Newtonian 
cooling  scheme  shows  a  5  degree  warming  at  the  top  of  the 
model  in  the  winter  polar  region.  By  equation  (11)  this 
appears  must  be  due  to  horizontal  advection  of  heat  since 
there  is  weak  radiative  and  adiabatic  vertical  motion 
cooling  in  this  region  (see  Figures  7  and  12)  throughout 
the  period .  At  40  km  at  the  winter  pole  there  has  been 
6  degree  warming  smoothing  the  temperature  minimum  found 
there  at  slightly  less  than  the  diabatic  heating  rate  of 
l°K/day.  There  is  a  major  20 °K  cooling  in  the  summer  trop¬ 
osphere  at  200  mb  in  both  integrations,  but  the  radiative 
cooling  rates  are  much  too  weak  to  have  played  a  major 
role  there . 

In  Figure  11  after  5  days  GCM  integration  using  the 
Ramanathan  model  there  is  little  change  in  temperature  in 
the  winter  polar  region  at  the  top  of  the  GCM.  In  compar¬ 
ison  with  the  Newtonian  model  run,  radiative  cooling  is 
greater  in  the  Ramanathan  model .  This  greater  cooling 


counteracts  the  horizontal  heat  convergence  noted  above. 

At  the  temperature  minimum  at  40  km  the  temperature  has 
slipped  four  degrees  to  an  even  colder  value.  In  both 
regions  the  band  model  radiational  cooling  has  increased 
the  anomalous  pole-to-pole  temperature  gradient  in 
comparison  with  the  Newtonian  scheme. 

One  effect  of  radiational  cooling  which  could  be 
expected  would  be  a  compensatory  adiabatic  wanning  by 
sinking  motions.  Figures  12  and  13  show  the  zonal  average 
vertical  velocities  field  found  in  the  model  after  5  days 
with  the  two  radiation  schemes.  There  are  no  significant 
differences  in  the  two  fields  suggesting  that  the  adiabatic 
warming  induced  by  the  increased  radiational  cooling  in 
the  winter  polar  region  is  negligible.  The  vertical  veloci¬ 
ties  are,  in  fact,  upward  suggesting  that  the  major  energy 
balance  is  between  horizontal  advection  of  heat  and  the 
adiabatic  cooling  by  the  vertical  velocities  which  are 
fairly  steady  over  the  5  day  period  above  30  km.  In  Tables 
2  and  3  are  the  magnitudes  of  the  terms  of  the  thermodynamic 
prediction  equation  for  the  upper  atmosphere.  Horizontal 
heat  advection  was  computed  as  a  residual  from  the  other 
components.  The  major  balance  terms  are  the  horizontal 
influx  of  heat  by  advection  and  the  adiabatic  cooling  by 
rising  vertical  motions.  Comparison  of  the  two  tables 
suggest  that,  at  least  in  the  initial  5  day  period  the  in¬ 
creased  cooling  by  the  band  model  over  the  Newtonian 
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TABLE  2 

Height 

70km 

60km 

50km 

40km 

30km 

25km 


:  Components  of  the  Thermodynamic  Heating  in  the 
Newtonian  Five  Day  Model  Run  at  80°N  (degK/day) 

3T  dZ 

Model  Level  “St  =  -k  x  Vf  *VT  ~dt 0 


2 

1.1 

12.6 

-10.9 

-  , 

3 

1.2 

21.3 

-18.6 

-1 

4 

1.0 

26.9 

-24.1 

-1 

5 

1.1 

31.6 

-28.5 

-2 

6 

1.1 

32.0 

-28.6 

-2 

7 

1.1 

31.0 

-27.4 

-2 

8 

1.1 

29.7 

-26.4 

-2 

9 

1.5 

33.2 

-29.5 

-2 

10 

1.6 

26.9 

-23.9 

-1 

11 

1.1 

20.3 

-19.4 

12 

1.4 

14.7 

-13.5 

13 

1.9 

13.0 

-10.4 

- 

14 

1.1 

10.5 

-  8.4 

-1 

15 

-.7 

7.0 

-  6.5 

-1 

16 

+  .8 

7.5 

-  4.7 

-2 

17 

.0 

5.4 

-  3.6 

-1 

trio 
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TABLE  3 . : 

Components 

of  the  Thermodynamic  Equation 

in  the 

Band  Model 

Five  Day 

GCM  Run. 

Height 

Model  Level 

3T 

St 

=  -k  x  VY  •  VT 

dZ 

~dt° 

70km 

2 

+  .2 

12.4 

-  9.7 

-2.5 

3 

-  .1 

20.2 

-16.8 

-3.5 

4 

-  .6 

26.7 

-23.5 

-3.8 

5 

-  .6 

31.9 

-27.7 

-4.8 

60km 

6 

-1.4 

33.1 

-28.2 

-6.3 

7 

-1.3 

32.2 

-27.4 

-6.1 

8 

-1.2 

30.9 

-26.7 

-5.4 

50km 

9 

-  .1 

34.5 

-30.3 

-4.3 

10 

-  .2 

28.2 

-24.6 

-3.8 

11 

-  .9 

22.1 

-20.3 

-2.7 

40km 

12 

-  .6 

15.2 

-14.1 

-1.7 

13 

+  .7 

12.2 

-10.4 

-1.1 

14 

-  .7 

9.4 

-  8.2 

-1.9 

15 

-1.3 

7.4 

-  6.7 

-2.0 

30km 

16 

+  .6 

7.2 

-  5.0 

-1.6 

25km 

17 

+  .3 

5.7 

-  3.6 

-1.8 

28 


scheme  primarily  affects  the  local  temperature  change.  The 
difference  in  the  local  temperature  changes  between  the 
model  runs  is  roughly  10%  of  the  horizontal  heat  advection 
term . 

There  is  some  dynamical  adjustment  to  the  increased 
radiational  cooling  even  in  the  five  day  model  run.  For 
example,  at  level  9  near  the  stratopause  the  horizontal 
and  vertical  adiabatic  terms  produced  4.2  deg  K/day  heating 
with  the  band  model.  This  was  .5  deg  K/day  stronger  than 
in  the  Newtonian  run.  This  reflects  an  increase  in 
circulation  intensity  with  a  negative  feedback  on  the 
increased  radiational  cooling  in  effect.  Longer  model 
runs  will  be  needed  to  determine  the  final  extent  of  the 
dynamical  adjustment  and  the  final  effects  on  the  meridional 
temperature  gradient. 

Manabe  and  Hunt  (1968,  Fig.  17)  determined  heat  bal¬ 
ance  components  in  an  18  level  stratospheric  general 
circulation  model  up  to  37.5  km.  At  79°N,  heat  was  added 
by  large  scale  eddies .  Heat  was  lost  by  the  mean  meridional 
circulation  (MMC)  and  radiative  effects  in  a  2:1  ratio  above 
25  km  in  their  model.  Table  3.  suggests  qualitative  agree¬ 
ment  with  Manabe  and  Hunt.  Vertical  velocities  are  weaker 
in  Manabe  and  Hunt's  model  in  the  25-40  km  range.  This 
is  probably  because  these  levels  are  near  the  rigid  upper 
lid  in  their  model. 

As  mentioned  in  the  discussion  of  the  thermodynamic 
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prediction  equation,  net  imbalances  in  the  horizontal  global 
heating  rates  are  effectively  subtracted  out  at  all  grid 
points.  These  net  imbalances  are  not  small  at  some  levels. 
Table  4  shows  the  net  global  heating  averages  for. the 
Ramanathan  model  heating  rates  in  Figure  8.  The  maximum 
value  occurs  at  level  10:  +  3.15  deg  K/day.  The  model 
equations  require  that  this  be  subtracted  out  before  entry 
into  the  prediction  equation.  So,  at  the  winter  pole, 
infrared  cooling  computed  as  +.31  became  -3.76  when  the 
global  heating  was  subtracted .  The  additional  error  of 
-,92°K  is  due  to  spectral  smoothing  of  the  UV  heating. 

Net  global  radiative  heating  in  the  upper  atmosphere 
has  been  noted  by  other  authors.  Johnson  and  Gottlieb 
(1971)  computed  vertical  eddy  heat  diffusion  coefficients 
as  a  function  of  altitude  based  on  the  assumption  that 
the  excess  global  heating  diffuses  to  lower  atmospheric 
levels  where  the  increased  densities  make  it  negligible. 

In  the  spectral  coordinates  a  simple  way  to  represent  this 
diffusion  would  be  to  adjust  the  UV  heating  components 
proportionally  to  balance  the  global  infrared  cooling  as 
shown  in  Table  4.  at  each  model  level: 


Q'  =  Q 

UV  UV 


QC02  +  QC>3 

Q 

uv 


(12) 


Net  global  heating  exceeds  1  deg/day  between  40  and 
60  km  as  shown  in  Table  4.  In  this  range,  the  adjustment 
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TABLE  5.: 


Band  Model  R.E.  Temperatures  and  GCM  Global 
Average  Temperatures. 


Model 

Level 

Approx . 

Height  (km) 

Radiative* 

Equilibrium 

Temperatures 

Global 

Average 

Temperatures 

1 

72 

•  242 

211 

3 

66 

248 

227 

5 

61 

263 

242 

7 

55 

276 

259 

9 

48 

282 

267 

11 

42 

260 

255 

13 

36 

238 

243 

15 

30 

226 

231 

17 

25 

220 

220 

19 

20 

216 

212 

21 

15 

220 

213 

23 

9 

236 

234 

25 

3 

276 

266 

*  Alimandi  and  Visconti  (1979) 
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factor  in  (12)  was  fairly  constant,  ranging  from  .65  -.79. 

An  analysis  of  net  diabatic  heating  rates  using  this 
adjustment  is  shown  in  Figure  9  .  In  comparison  with 
Figure  8,  the  unadjusted  band  model  heating  rates,  the 
general  shape  of  the  heating  analysis  appears  to  be  pre¬ 
served  but  the  magnitude  of  the  gradients  has  been  reduced. 

Net  global  heating  was  found  at  essentially  all 
model  levels  in  the  upper  atmosphere.  The  radiational  band 
model  is  trying  to  heat  the  entire  atmosphere  in  this 
region.  Alimandi  and  Visconti  (1979)  have  calculated  a 
radiative  equilibrium  temperature  profile  for  their  version 
of  the  band  model .  Since  the  band  model  reported  here  is 
based  on  theirs,  their  equilibrium  profile  should  apply  at 
least  approximately.  Table  5  showed  that  at  almost  all 
levels,  the  radiative  equilibrium  temperature  is  warmer 
than  the  specified  horizontal  average  temperatures  of  the 
GCM.  The  GCM  horizontal  average  temperatures  are  based  on 
observations,  but  the  upper  atmosphere  observations  are 
not  considered  highly  reliable.  A  straight  forward  method 
of  tuning  the  models  into  agreement  would  be  to  use  the 
radiative  equilibrium  temperatures  as  the  specified  hori¬ 
zontal  average  temperatures  in  the  GCM. 

6.  SUMMARY 

The  revised  Ramanthan  band  model  developed  in  this 
work  appears  to  be  in  good  agreement  with  the  work  of 


Dopplick  (1971),  Slade  (1969),  and  Park  and  London  (1974), 
while  being  much  faster  than  the  original  code. 

When  incorporated  into  the  MIT  stratospheric  GCM,  the 
band  model  created  much  stronger  pole-to-pole  differential 
heating  than  the  Newtonian  heating  scheme.  After  only  a 
five  day  period  of  integration,  the  band  model  sharpened 
the  existing  erroneous  temperature  gradients  by  several 
degrees.  However  it  was  noted  that  the  differences  in  local 
temperature  changes  between  the  heating  schemes  are  only 
10%  of  the  horizontal  heat  advection.  Poleward  horizontal 
eddy  heat  fluxes  tend  to  increase  with  increasing  resolu¬ 
tion  in  GCMs .  Thus  when  the  MIT  stratospheric  model  is 
converted  from  six  to  eighteen  wave  resolution,  the  band 
model  may  handle  heating  more  accurately  than  the  Newtonian 
model .  Dynamical  time  constants  are  longer  than  five  days 
and  some  dynamical  adjustment  was  noted.  Longer  model  runs 
will  be  appropriate  in  the  eighteen  wave  GCM  to  determine 
steady  state  average  heating  rates  and  meridional  tempera¬ 
ture  gradients . 

The  existence  of  large  net  global  heating  rates, 
especially  from  40-60  km,  was  identified  as  a  problem  in 
applying  the  radiation  band  model  to  the  balance  equations 
of  the  MIT  stratospheric  GCM.  To  adjust  the  UV  heating 
alone  would  require  a  20-35%  reduction  in  UV  heating  in 
the  40-60  km  region.  A  more  general  procedure  to  bring  the 
band  model  into  balance  in  the  GCM  would  be  to  use  the 
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band  model  radiative  equilibrium  temperatures  to  specify 
the  horizontal  average  temperatures  of  the  GCM. 
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Figure  7  s  Zonal  averaged  net  diabatic  heating  rates  with 
Newtonian  cooling,  Day  361,  Run  34. 
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Zonal  averaged  net  diabatic  heating  rates  with 
band  model  cooling,  Day  361,  Run  34. 
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Figure  9s  Zonal  averaged  net  diabatic  heating  rates  with 
adjusted  cooling,  Day  361,  Run  34. 
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averaged  temperature  field  after  5  days 
integration  with  band  model  cooling. 
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